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ABSTRACT 

Recent observations of high redshift quasar spectra reveal long gaps with little flux. A small 
or no detectable flux does not by itself imply the intergalactic medium (IGM) is neutral. Inferring 
the average neutral fraction from the observed absorption requires assumptions about clustering 
of the IGM, which the gravitational instability model supplies. Our most stringent constraint on 
the neutral fraction at z ~ 6 is derived from the mean Lyman-beta transmission measured from 
the z = 6.28 SDSS quasar of Becker et al. - the neutral hydrogen fraction at mean density has to 
be larger than 4.7 x 10^"*. This is substantially higher than the neutral fraction of ^ 3 — 5 x 10~^ 
at 2; = 4.5 — 5.7, suggesting that dramatic changes take place around or just before z ~ 6, 
even though current constraints are still consistent with a fairly ionized IGM at 2 ~ 6. These 
constraints translate also into constraints on the ionizing background, subject to uncertainties 
in the IGM temperature. An interesting alternative method to constrain the neutral fraction 
is to consider the probability of having many consecutive pixels with little flux, which is small 
unless the neutral fraction is high. It turns out that this constraint is slightly weaker than the 
one obtained from the mean transmission. We show that while the derived neutral fraction at 
a given redshift is sensitive to the power spectrum normalization, the size of the jump around 
2; ~ 6 is not. We caution that the main systematic uncertainties include spatial fluctuations 
in the ionizing background, and the continuum placement. Tests are proposed. In particular, 
the sightline to sightline dispersion in mean transmission might provide a useful diagnostic. We 
express the dispersion in terms of the transmission power spectrum, and develop a method to 
calculate the dispersion for spectra that are longer than the typical simulation box. 



Subject headings: cosmology: theory - intergalactic medium - large scale structure of universe; 
quasars - absorption lines 



1. Introduction 

Recent spectroscopic observations of z <; 4.5 quasars discovered by the Sloan Digital Sky Survey (SDSS) 
have opened up new windows into the study of the high redshift intergalactic medium (IGM) (Fan et al. 
2000, Zheng et al. 2000, Schneider et al. 2001, Anderson et al. 2001 Fan et al. 2001a, Becker et al. 2001, 
Djorgovski et al. 2001). In particular, Becker et al. (2001) observed Gunn-Peterson troughs (Gunn & 
Peterson 1965) in the spectrum of a z = 6.28 quasar, which were interpreted as suggesting that the universe 
was close to the reionization epoch at z ~ 6. 
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That the absorption increases quickly with rcdshift is not by itself surprising: ionization equilibrium 
tells us that the neutral hydrogen density is proportional to the gas density squared, which is proportional to 
(1 + z)^ at the cosmic mean. The evolution of the ionizing background and gas temperature will modify this 
redshift dependence, but the rapid evolution of absorption remains a robust outcome. What is interesting, as 
Becker et al. (2001) emphasized, is that the observed mean transmission at redshift z '--^ 6 is lower than what 
one would expect based on an extrapolation of the column density distribution and its redshift evolution 
(number density of clouds scaling as ~ (1 + z)^'^) from lower redshifts. On the other hand, the popular 
gravitational instability theory of structure formation provides detailed predictions for how the IGM should 
be clustered, and how this clustering evolves with redshift, which has been shown to be quite successful 
when compared with observations at z ~ 2 — 4 (See e.g. Cen et al. 1994, Zhang et al. 1995, Reisenegger 
& Miralda-Escude 1995, Hernquist et al. 1996, Miralda-Escude 1996, Muecket et al. 1996, Bi & Davidsen 
1997, Bond & Wadsley 1997, Hui et al. 1997, Croft et al. 1998, Theuns et al. 1999, Bryan et al. 1999, 
McDonald et al. 2000a). These predictions allow us to directly infer the neutral fraction of the IGM from 
the observed absorption (the relation between the two depends on the nature of clustering of the IGM), and 
so can further inform our interpretations of the recent z 6 results. 

How neutral is the IGM at z ~ 6, and how different is the neutral fraction compared to lower redshifts? 
These are the questions we would like to address quantitatively, making use of the gravitational instability 
model of the IGM. 

The paper is organized as follows. First, we start with a brief description of the gravitational instability 
model for the IGM and the simulation technique in §|^. In § |3.1^ we derive the neutral hydrogen fraction 
Xhi, and equivalently the level of ionizing flux Jhi, at several different redshifts leading up to z 6 from 
the observed mean Lyman-alpha (Lya) transmission. This exercise using the Lya spectrum is similar to the 
one carried out in McDonald & Miralda-Escude (2001), except for the addition of new high redshift data. 
Q We then examine in §|3.2| the constraints on the same quantities Xm and Jhi from the observed mean 



Lyman-beta (Ly/3) transmission, Ly/3 being particularly useful at high Lya optical depth, because the Ly/3 
absorption cross-section is a factor of ~ 5 smaller than the Lya cross-section. The goal here is to use Ly/3 



absorption to obtain constraints on X^n and Jhi that are as stringent as possible. In §3.2, we also examine 
the sensitivity of our conclusions to the power spectrum normalization. 

An intriguing question is: instead of focusing on the mean transmission, can one make use of the fact 
that the observed spectrum at z '--^ 6 contains a continuous and long stretch (^ 200 — 300A) with little or no 
detected flux to obtain more stringent limits on the neutral fraction or Jhi? The idea is that since the IGM 
gas density naturally fluctuates spatially, it seems a priori unlikely to have no signiflcant upward fluctuation 
in transmission for many pixels in a row - unless of course the neutral fraction Xm is indeed quite high. 



We will show in §3.3 this provides constraints that are slightly weaker to those obtained using the mean 
transmission. 

In all the simulations discussed in this paper, the ionizing background is assumed uniform spatially, 
just as in the majority of high redshift IGM simulations. A natural worry is that as the universe becomes 
more neutral at higher redshifts, the ionizing background would be more non-uniform. One way to test this 
is to use several lines of sight, available at z ~ 5.5, and compare the observed line-of-sight scatter in mean 



^This part of the calculation involving the matching of the mean Lya transmission is also similar to a number of earlier 
papers where the primary quantity of interest is the baryon density (e.g. Ranch et al. 1997, Weinberg et al. 1997, Choudhury 
et al. 2000, Hui et al. 2001). Here, we fix the baryon density and study the ionizing background or the neutral fraction instead 
(see §|). 
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transmission against the predicted scatter based on simulations with a uniform background. We discuss this 
in estimate the level of ionizing background fluctuations, and make predictions for the scatter at z ^ 6. 
Here, we also introduce a technique to handle the problem of limited box-size. 

Readers who are not interested in details can skip to §H where we summarize the constraints obtained. 
We also discuss in §|| the issue of continuum placement, and how the associated uncertainties can be es- 
timated. While the work described in this paper was being carried out, several papers appeared which 
investigate related issues (Barkana 2001, Razoumov et al. 2001, Cen & McDonald 2001, Gnedin 2001, Fan 
et al. 2001b). Where there is overlap, our results are in broad agreement with these papers. We present 
a comparision with other authors at the end of §3.S. Our approach here is most similar to that of Cen & 
McDonald (2001). In addition to obtaining constraints on the ionizing background from the Lya and Ly/3 
transmission as was considered by Cen & McDonald, we consider the possible constraint from the Gunn- 
Peterson trough itself, examine the dependence on power spectrum normalization, and develop a method to 
predict the scatter in mean transmission by relating it to the power spectrum, which might be of wider in- 
terest. We also place a slightly stronger emphasis on the neutral fraction, which is more robustly determined 
compared to the ionizing background or photoionization rate. 



2. The Gravitational Instability Model for the IGM 

The Lya optical depth is related to the IGM density, assuming ionization equilibrium , via 

r„=A„(l + 5)2-0-7(7-1) (1) 

where S is the gas overdensity [6 ~ {p ~ p)/P7 where p is the gas density and p its mean), 7 is the equation 
of state index for the IGM^, and Aa is given by (see e.g. Hui et al. 2001 and references therein): 
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where Xm = nm/nn (rifj is the total density of neutral and ionized hydrogen, and nm is the neutral 
hydrogen density) is the neutral hydrogen fraction at mean density ((5 = 0). ^ Here H{z) is the Hubble 
parameter at redshift z, Hq is the Hubble parameter today, Hq = IQQh km/s/Mpc, fib is the baryon density 
in units of the critical density. The value of 11.7 for H{z)/Hq above corresponds to that appropriate for a 
cosmology with fi^ = 0.4 and JIa = 0.6 at z = 6, where flm and JIa are the matter and vacuum densities 
in units of the critical density today. 



^ The photoionized IGM at overdensity of a few or less is expected to follow a tight temperature-density relation of the form 
T = To(l + (5)7— 1 , where T is the gas temperature and Tq is its value at the cosmic mean density (see Hui & Gnedin 1997). We 
caution that close to reionization, these quantities may not be a function of 5 alone. The IGM may bo heated inhomogenously, 
causing spatial fluctuations in To and 7. 

^The neutral fraction at arbitrary 5 is given by Xhi times (1 -I- (5)i-0-^('''-l) Throughout this paper, whenever we quote 
values for Xjji, we refer to the neutral hydrogen fraction at the cosmic mean density (5 = 0. 
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The neutral fraction Xm is related to the ionizing background by ^ 



Xm = 1.6 X 10" 
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where the dimensionless quantity Jjji is related to the photoionization rate Ffji by 

Thi = 4.3 X lO^^Vms^i 



(3) 



(4) 



The quantity Jhi provides a convenient way of describing the normalization of the ionizing background, 
without specifying the exact spectrum, in a way that is directly related to the physically relevant Fhi (e.g. 
Miralda-Escude et al. 1996). It is related to the specific intensity at 912A ji,^^ by Jhi = Ji/hi x [3/(/3 + 3)], 
where /3 is the slope of the specific intensity just blueward of 912 A {ju oc where v is frequency), and ju^^i 
is measured in the customary units of 10~^^ ergs^^ cm~^ Hz~^ ster~^ (for non-power law j^, eq. provides 
the exact definition for Jhi; see e.g. Hui et al. 2001). 



Two more ingredients should be mentioned to complete the specification of our model for the Lya 
absorption (see e.g. Hui et al. 1997 for details). First, the optical depth as a function of velocity is computed 
by taking the right hand side of eq. (^ in velocity space (i.e. taking into account peculiar velocities) and 
smoothing it with a thermal broadening window. Second, the gas density and velocity fields are predicted 
by some Cold Dark Matter (CDM) cosmological model using numerical simulations. 

There are obviously a number of free parameters in our model. Let us discuss each of them in turn. 

Throughout this paper, we assume fib/i^ = 0.02, as supported by recent cosmic microwave background 
measurements (Netterfield et al. 2001, Pryke et al. 2001) and the nucleosynthesis constraint from primordial 
deuterium abundance (Buries et al. 2001). We also assume throughout h — 0.65, Vim — 0.4, and JIa = 0.6. 
Variations of these parameters within the current bounds do not contribute significantly to the uncertainties 
of the constraints obtained in this paper (see Hui et al. 2001). 

The temperature Tq and equation of state index 7 at the redshifts of interest in this paper are somewhat 
uncertain. There are no direct measurements of the thermal state of the IGM at our redshifts of interest, 
2; ^ 4. Measurements at 2: ^ 4 yield values consistent with Tq = 2 x 10^ K and 7 = 1 (McDonald et al. 
2000b, Ricotti et al. 2000, Zaldarriaga et al. 2001. Schaye et al. 2000, however, measure a slightly lower 
temperature). Given that the temperature right after reionization is expected to be about 25000 K with 
7 = 1 (with some dependence on the hardness of the ionizing spectrum; see e.g. Hui & Gnedin 1997), which 
is not too different from the measurements at z ^ 4, we will assume throughout this paper, when making 
use of eq. (||) to infer Jhi, that Tq = 2 x 10* K and 7 = 1. Note that while the theoretically allowed range 
for 7 is from 1 to 1.6 (Hui & Gnedin 1997), what matters for our purpose is 2 — 0.7(7 ~ 1) (sq. lll|)i which 
only ranges from 1.58 to 2, and does not significantly affect our results. It is also important to emphasize 
that the inference 0/ Xhi from observations, unlike the case for Jhi, is not subject directly to uncertainties 
in the temperature Tq. This is because observations constrain Aa from which we can obtain Xm without 
knowing Tq (see eq. [D). 



^This equation assumes that hydrogen is highly ionized and that hchum is largely doubly ionized. If helium is only singly 
ionized, the relation between Jjji and Xjji will be changed slightly: the right hand side of eq (js]) will be multiplied by 0.93. 

^ The above statement is subject to two small caveats. First, the optical depth given in eq. (|l|) has to be smoothed with a 
thermal broadening window whose width depends on Tq. We find that in practice, the exact width of the thermal broadening 
kernel does not affect very much quantities such as the mean transmission, which is what we will be interested in. Second, Tq 
also affects the gas dynamics via the pressure term in the equation of motion. As we will discuss below, the effect of varying 
To also appears to be small in this regard. 
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To generate realizations of the density and velocity fields for a given cosmology, we run Hydro-Particle- 
Mesh (HPM) simulations (Gnedin & Hui 1998). The HPM algorithm is essentially a Particle-Mesh code, 
modified to incorporate a force term due to gas pressure in the equation of motion. ^ For the initial power 
spectrum, we use a Cold Dark Matter (CDM) type transfer function, as parameterized by Ma (1996), which 
is very similar to the commonly used Bardeen et al. (1986) transfer function. For the primordial spectral 
slope, we adopt n = 0.93 (Croft et al. 2000, McDonald et al. 2000a). For the linear power spectrum 
normalization, we employ the range suggested by measurements from the Lya forest of Croft et al. (2000): 
A2(fc) EE 47rfc3p(fc)/(27r)3 = O.lAtoil at z = 2.72 at a velocity scale of fc = O.G3( km/s)-i.0 We, however, 
caution that the error-bar given is somewhat dependent on the assumed error of the mean transmission 
measurements, which is sensitive to the accuracy of the continuum-fitting procedure (see e.g. Zaldarriaga 
et al. 2001 for a slightly different assessment of the error-bar). The power spectrum in this model has a 
similar shape to that of favored cosmological models, but slightly lower amplitude (Croft et al. 2000). In 



§3.2 we demonstrate that our main conclusion, that the neutral fraction increases dramatically near z ^ 6, is 
insensitive to our assumptions about the amplitude of the power spectrum. In practice, we examine models 
with different normalizations by running a simulation with outputs at several different redshifts: each redshift 
then corresponds to a different power spectrum normalization, and linear interpolation is performed to reach 
any desired normalization. ^ 

Our simulations have a box size of 8.9 Mpc/h, with a 256'^ grid. McDonald & Miralda-Escude (2001) 
found this box size and resolution to be adequate for IGM studies up to z ~ 5. We have verified that the 
same is true up to z = 6, in the sense that the transmission probability distribution has converged for our 
choice of simulation size and resolution. 

Finally, we should say a few words about simulations of the Ly(3 region. In regions of the quasar 
spectrum that are between 973A(1 -|- Zem.) and 1026A(1 -I- Zem.), where Zem. is the redshift of the quasar, 
two kinds of absorption can exist: one is Ly/3 due to material at redshift 0.948(1 -I- Zom.) < 1 -I- z < 1 -I- Zem., 
the other is Lya due to material at redshift 0.800(1 -I- Zom.) < 1 + z < 0.844(1 -I- Zcm.). In other words, in 
such a region, the observed optical depth would be given by t = Tfj + where Tp and Tq, arises at different 
redshifts. The Lya optical depth can be computed as before. The Ly/3 optical depth Tp can be computed 
using eq. (|l|), except that is replaced by Ap: 

A, = (5) 

The factor of 5.27 reflects the fact that the Ly/3 transition has a cross-section that is 5.27 times smaller than 
Lya. 



®The temperature-density relation has to be specified as a function of redshift in the HPM code to compute the pressure 
term. We follow McDonald & Miralda-Escude (2001) and linearly interpolate between Tq = 19000 K and 7 = 1.2 at 2 = 3.9 
and To = 25000 K and 7 = 1 at the redshift of reionization Zrcion- We found that assuming Zreion = 7 versus Zj-^i^n = 10 
results in negligible difference in our results, in particular concerning the mean decrement and the probability distribution of 
transmission. All results in this paper are quoted from the Zreion = 7 HPM simulations. Note that in inferring Xjii and Jjii 
from eq. and we always use Tq = 20000 K and 7 = 1 for simplicity, as mentioned before. 

'^This normalization corrects an error in an earlier draft of Croft ct al. (2000). (R. Croft, private communication.) 

®We do not vary the primordial spectral index n here. Quantities such as the mean transmission which we are interested in 
here are generally sensitive to power on only a small range of scales. Varying n is therefore largely degenerate with varying A^. 
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3. Constraints on the Neutral Hydrogen Fraction and the Ionizing Background 

3.1. Constraints from the Lya Mean Transmission 

Using eq. and (||), we compute the Xm, which also fixes Jhi (eq. [||), necessary to match the 
observed Lya mean transmission (e~'^°) at z — 4.5 — 6 (see Table 1 for a summary of the measurements). 
The results of our calculation are presented in Fig. 0. This plot also contains a point a,t z — 6.05 which is 
the result of matching the mean transmission in the Ly/3 forest, as we describe in § ^.2| . 

Also shown in the figure is a dotted line which shows X^ii oc {1 + zY, which appears to be a good fit 
to the data from z — 4.5 to z = 5.7. From eq. (||), one can see that such a trend for the neutral fraction is 
equivalent to assuming constant Jhi (or more accurately, constant JhiTq '^; see eq. [^). 

As one can see, ignoring for now the Ly/? point, the neutral fraction does appear to have a modest jump 
around z 6: it increases by a factor of 4.0 from z ~ 5.7 to z = 6.05, while it changes at most by ~ 1.9 
from z = 4.5 to z = 5.7. A similar trend (but opposite in sign) can be seen in the ionizing flux Jhi- The 
1 a error-bar here takes into account the measurement error in mean transmission, and the range of power 
spectrum normalization stated in As we have explained in while Xjji is not sensitive to the assumed 
temperature of the IGM (Tq); our constraints on Jhi are directly influenced by it. As emphasized before, 
we assume Tq = 2 x 10* K. In other words, our constraints on Jhi are really constraints on the quantity 
>^Hi(To/2 X 10*K)°''^ (see eq. [||). It is therefore straightforward to rescale our constraints on Jhi if the 
temperature were a little bit different. ^ It is an interesting question to ask whether the apparent jump in 
the ionizing flux can instead be attributed to a jump in the temperature. In general, the temperature is 
expected to evolve slowly with redshift after reionization (Hui & Gnedin 1997). 

Regarding the measurement error, we should also emphasize that the Becker et al.'s 2 a error actually 
includes the possibility of having zero transmission at z = 6.05. This means that at 2 cr, we would only 
have a lower limit on Xhi, or an upper limit on Jhi, for the highest redshift point in Fig. |l|, allowing the 
possibility that the IGM is neutral at z ^ 6, Xhi — 1. 



3.2. Constraints from the Ly/3 Mean Transmission 

In this section we consider the constraints placed by Becker et al.'s measurement of the mean trans- 
mission in the Ly/3 region. Absorption in the Ly/3 region has two components: t — t^A- t^, where the Lya 
optical Tq and the Ly/3 optical depth originate at different redshifts. Ly/3 absorption due to material at 
z = 6 coincides in wavelength with Lya absorption due to material at z = (1 -f- 6) x 1026/1216 — 1 = 4.9. Be- 
cause the points of origin are so widely separated, they can be effectively treated as statistically independent 
i.e. (e^'^) — (e~'^°)(e~'^*'). Becker et al. measured (e""^*^) at z ~ 6 by dividing the net mean transmission 
[eT'^) in the Ly/3 region by the mean transmission in Lya at z ^ 5. They obtained (e^'^f ) = —0.002 ± 0.020. 
Clearly, this measurement is consistent with a completely neutral IGM. However, the interesting question is: 
what kind of lower limit does it set on the neutral fraction, and does it improve upon the lower limit from 
the mean Lya absorption ? 

We carry out a calculation that is analogous to what is described in § |3.lt except for the key difference 



^The temperature also affects the thermal broadening window, but wo find that in practice its effect on our constraints on 
Aa (eq. [pj) is small; see discussion in 33. 
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that in computing from the simulated density and velocity fields, we employ Ap which is a factor of 5.27 
smaller than Aa (see eq. & [||). The results of our calculation are shown in Fig. |l|, as the highest redshift 
points in the plot, which have error bar arrows pointing towards a completely neutral IGM and a vanishing 
ionizing background. It can be seen that the (1 — cr) lower limit on Xhi is Xm > 4.7 x 10^''. This is a factor 
of ^ 3 larger than the neutral fraction required to match the upper limit on the mean transmission in Lya 



for our fiducial normalization, and a slightly stronger constraint than that obtained in section 3.1, including 
the uncertainty in power spectrum normalization. Similarly, the upper limit on Jhi is Jhi < 9.0 x 10^^. 
The moral here is that because the Ly/3 absorption cross-section is a factor of 5.27 smaller than the Lya 
cross-section, Ly/3 offers a more sensitive probe of the neutral fraction, especially when the Lya optical depth 
is high. 

The neutral fraction at z ~ 6 is thus a factor of ^ 10 higher than that at redshift z ^ 5.7, where it is 
Xm = 4.9 X 10"^. This dramatic change in the neutral fraction is suggestive, probably indicating that the 
rcionization epoch is nearby. 

Furthermore, this conclusion is not sensitive to our assumptions about the amplitude of the power 
spectrum. Although the neutral fraction at redshift z = 6.05 is itself sensitive to the amplitude of the power 
spectrum, we find that the factor by which the neutral fraction increases from z = 5.7 to z — 6.05 depends 
only weakly on the amplitude. In Fig. || we plot both the neutral fraction at z = 6.05 and the jump in the 
neutral fraction for a range of different power spectrum normalizations. The jump is defined as the ratio 
Xm{z — 6.05)/Xhi{z — 5.7). Here Xhi{z = 6.05) is the lower limit resulting from the 1 a error in the 
mean transmission in Ly/3 at z = 6.05 and the error bars in the jump arise from the 1 a error in the mean 
transmission in Lya at z = 5.7. As one can see in the plot, the lower limit on the neutral fraction at z — 6.05 
varies from Xm > 3 x 10"* to Xm > 9 x lO"'^ as A'^{k = 0.03(km/s)~\ z = 2.72) varies from 0.5 to 1.3. 
The neutral fraction itself varies significantly with power spectrum normalization, scaling approximately as 
Xm oc [A^(fc — 0.03s/km, z = 2.72)]^-^, for this range of normalizations. The jump, however, changes only 
slightly over a large range of normalizations. As A^(fc — 0.03( km/s)"^, z — 2.72) varies from 0.5 to 1.3, 
the jump changes only from ~ 9.7 to ^ 11.1. Our conclusion that the neutral fraction of the IGM increases 
dramatically near z 6 seems robust. 

One can also consider the absorption in the Ly7 region, or even the higher Lyman series. In practice, 
the accumulated amount of absorption from Lya as well as Ly/3 at different redshifts makes it harder to 
measure the Ly7 transmission itself with good accuracy. 

Our constraints on the neutral fraction and the intensity of the ionizing background are consistent 
with those found by other authors, given our different choices of power spectrum normalization. Fan ct 
al. (2001b) found, from the mean Ly/3 transmission, that r_i2 < 0.025, where r_i2 is the photoionization 
rate of eq (^) in units of 10~^^s~^. Although this constraint is somewhat stronger than the constraint 
implied by our fiducial model, r_i2 < 0.039, we expect the difference due to our different power spectrum 
normalizations. Fan et al.'s (2001b) constraint comes from semi-analytic arguments, shown consistent with 
an LCDM simulation with n„, = 0.3, flA = 0.7,/i = 0.65, Qbh^ = 0.02, and as = 0.9. This model has 
a substantially larger normalization, A^(fc = 0.03( km/s)~^, z — 2.72) — 1.25, than our fiducial model of 
A^(fc = 0.03( km/s)"^, z ~ 2.72) = 0.74. The difference in normalization refiects some tension between the 
normalization derived from the observed cluster abundance, which Fan uses, and that from the Lyman-a 
forest which our model is based on (Croft et al. 2000). Fan et al. (2001b) assume Tq = 2.0 x 10* K in 
placing their constraint. Their limit, r_i2 < 0.025, includes only uncertainties in the mean transmission 
and not additional uncertainties from the power spectrum normalization. From figure (H), we infer that Fan 
et al.'s (2001b) normalization implies X^ii > 8.8 x 10"* in our cosmology. Rescaling this result from our 
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assumed ftm — 0.4 to an ~ 0.3 cosmology, and using eqs. (^) and (^) we predict r_i2 < 0.024, or 
Xhi > 7.6 X 10^"*, for Fan et al.'s (2001b) model. The constraint of Fan et al. (2001b) is thus consistent with 
our constraint given our different choice of normalization. Cen & McDonald (2001), using a model similar 
to that of Fan et al. (2001b), obtained the constraint r_i2 < 0.032, using the Ly/3 mean transmission. The 
constraint is slightly weaker than that of Fan et al. (2001b) and our extrapolation to their normalization, 
because Cen & McDonald (2001) consider a larger upper limit to the observed mean transmission, including 
an estimate of the uncertainty due to sky subtraction. At slightly lower redshifts, we can also compare 
with the results of McDonald & Miralda-Escude (2001) derived from matching the mean Lya transmission. 
For example, aX z — 5.2, these authors found r_i2 = 0.16 to match the observed mean transmission of 
(e"'^'') = 0.09. McDonald & Miralda-Escude (2001) consider a model whose normalization we infer to be 
A^(fc = 0.03( km/s)^^, 2 — 2.72) — 0.98. To match the same mean transmission with this normalization 
we infer a somewhat higher photoionization rate, r_i2 = 0.19. Part of the difference may be that the 
r_i2 necessary to match a given mean transmission varies by ~ 5% between two different realizations of 
the density field. The remaining difference may come from the procedure of linearly interpolating between 
outputs or from some modeling difference. At any rate, our results are roughly consistent with those of other 
authors given our different power spectrum normalizations. 



3.3. Constraints from the Gunn-Peterson Trough Itself the Fluctuation Method 

The fact that Becker et al. (2001) observed a Gunn-Peterson trough, where a long stretch of the 
spectrum contains little or no flux, can conceivably be used to further tighten the constraints obtained from 
the previous sections. Since the IGM is expected to have spatial fluctuations, the probability of having many 
pixels in a row turning up a very small transmission must be low, unless the neutral fraction is intrinsically 
quite high. The same reasoning can be applied to either the Lya or Ly/3 absorption. We will discuss our 
method for Lya in detail. The method for Ly/3 is a straightforward extension. For simplicity, we will call 
this method, the fluctuation method. 

Becker et al. (2001) finds from the spectrum of SDSS 1030+0524, the z = 6.28 quasar, the Lya 
transmission is consistently below about 0.06 for a region that spans 260A, between 8450A to 8710A. The 
noise level per 4A pixel is -^Z (n^) ~ 0.02, where n represents the photon noise fluctuation. ^ The observed 
transmission F at a given pixel is F = e^'^ + n, where e^"^ is the true transmission. The noise here should be 
dominated by Poisson fluctuations of the subtracted sky background (as well as perhaps read-out error). Let 
P{Fi,F2, ...Fis[)dFi...dFN be the probability that N consecutive pixels have observed transmission fall into 
the range Fi ± dFi/2 ... F^ ± dFi^ /2. In our case, iV = 65 for the pixel size of AA. The problem is then to 
find the probability gg ... P{Fi...Fj^)dFi...dFN as a function of Jhi, and ask what maximal Jhi (or 

equivalently, minimal Xri) would give an acceptable probability. By choosing the "acceptable probability" 
to be within 68% of the maximum likelihood (maximum likelihood is achieved when the neutral fraction is 
unity), we obtain the 1 a upper limit on Jhi or 1 a lower limit on A"hi. 

Our simulation has a comoving box size of 8.9 Mpc/h, corresponding to A2A for Lya at z 6, which 
falls short of the wavelength range we need for this problem, which is 260A. In other words, the probability 



'^"We estimate the noise per pixel from Becker et al.'s error-bar in the mean transmission, which is ~ 0.003. This is 
estimated from a chunk of the spectrum which is 260A long, and so the dispersion per 4A pixel should be approximately 
{rfi) ~ 0.003 X v^65 ~ 0.02. Note that the actual dispersion varies across the spectrum, but this should suffice as a rough 
estimate. This estimate also agrees with an estimate of the error by comparing Fig. 1 and Fig. 3 of Becker et al. 
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P{Fi, F2, ■■■FN)dFi...dFN can be estimated directly from our simulation only for N < 10. However, the 
mass correlation length scale at this redshift ( <^ 1 Mpc/h) is actually a fraction of the box size, which means 
one can treat fluctuations on scales beyond the box size as roughly uncorrelated. Assuming so, we estimate 
/<o 06 ■■■ /<o 06 ^^^1' ■•■^io)'^-P'i---<^^io from the simulation, and then take its 6-th power, which would 
give us the probability that 60 consecutive pixels have transmission below 0.06. This is slightly smaller than 
the number 65 that we need, but at least will provide us conservative constraints on Xm and Jhi. We have 
also tested our approach by using fractions of the box-size as a unit, and find that our results do not change 
significantly (less than 10%). 

Fig. ||a (dotted curve) shows our estimate of the probability J^^ •■• /<o 06 ^(-^i' -^2, ...-FV)d-Fi...dFjv 
for = 60 and pixel size 4A, as a function of Jhi. Our simulated spectra have been convolved with 
the observation resolution (fuU-width-at-half-maximum of 1.8 A), rebinned into pixels of iA each and added 
Gaussian noise with a dispersion of 0.02. From the dotted curve in Fig. applying a likelihood analysis, we 
obtain a 1 ct upper limit on Jhi of Jhi < 0.014, and a corresponding lower limit on Xm of Xhi > 2.95 x 10^**. 
This is for a model with a power spectrum normalization of A^(fc = 0.03s/km, z = 2.72) — 0.74 (see 
The mean Lya transmission constraints for the same model are Jhi < 0.028 and Xhi > 1.5 x lO^"'. ^ This 
means that considering Lya alone, the fluctuation method yields somewhat stronger constraints compared to 
using simply the mean transmission. 

Fig. I^b (dotted curve) shows the same methodology applied to the Ly/3 Gunn-Peterson trough. A new 
ingredient here is that one needs an additional simulation of the same model at redshift z — 4.9 to produce the 
Lya absorption that can be overlaid on top of the Ly/3 absorption from z = 6.05. This additional simulation 
should have different initial phases to mimic the fact that fluctuations aX z ~ 4.9 and those a.t z — 6.05 
should be uncorrelated. We obtain 1 a limits of Jhi < 0.012 and Ahi > 3.4 x 10^*. This is about 40% weaker 
than the constraints we obtain from the Ly/3 mean transmission. In other words, from Lyj3 absorption, the 
fluctuation method yields slightly weaker constraints compared to using the mean transmission. It is also only 
slightly stronger than the constraint obtained from the fluctuation method applied to Lya. 

It is an interesting question to ask how many sightlines one would need to improve the constraints by, 
say a factor of 2. Our approach can be easily extended to multiple (uncorrelated) sightlines, and we find that 
about 5 sightlines (each containing a Gunn-Pcterson trough of the same length and same signal to noise) are 
necessary for such an improvement. 

Part of the difficulty with obtaining stronger constraints, in addition to the small number of sightlines, 
is the dominance of noise. The lower panel of Fig. ^ shows the one-pixel (4A) probability distribution 
function (PDF) of the true transmission e^'^ (i.e. no noise added), for three different values of Jhi (the 
power spectrum normalization is the same as that in Fig. The upper panel shows the corresponding 
PDF's of the observed transmission F (i.e. after convolving P{e^'^) with a Gaussian of dispersion 0.02). As 
expected, noisy data make the PDF's more similar. Nonetheless, as we pointed out above, with sufficient 
number of sightlines, there might be a non-negligible chance of seeing pixels with high transmission that take 
place at the tail of the PDF's, hence allowing us to distinguish between the different levels of the ionizing 
background. Alternatively, one can try improving the signal-to-noise per pixel. In Fig. ^jb, we show with a 
dashed curve the corresponding probability if the noise per pixel is lowered by a factor of 4. The constraints 
improve by a little more than a factor of 2. We should emphasize, however, systematic errors are likely 



^^Do not confuse these constraints, which are for the particular power spectrum normahzation mentioned above, to the 
constraints discussed in earUer sections, which include the uncertainty in the power spectrum normalization. We focus on a 
single model in this section for simplicity. 
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important here - we will discuss them in the next two sections. 

4. The Variance of the Mean Transmission 

If, as is suggested by our discussion in j ^3.2| (see Fig. |^), the IGM is close to the epoch of reionization 
at z ^ 6, one might expect large fluctuations in the ionizing background near that time. For instance, one 
line of sight might probe a region of the IGM where the ionized bubbles around galaxies or quasars have 
percolated, while another might probe the pre-percolation IGM. As mentioned before, the simulations we 
employ do not take into account fluctuations in the ionizing background. (For simulations incorporating 
radiative transfer see e.g., Gnedin & Abel 2001, Razoumov 2001). One useful check would then be to predict 
the sightline to sightline scatter in mean transmission from our simulations, and compare that against the 
observed scatter. At z ~ 5.5, 4 lines of sight are available for a measurement of the scatter. We will examine 
this, as well as make predictions for the scatter at z ~ 6, which more high redshift quasars in the future will 
allow us to measure. 

Our estimate relies on simulation measurements of the transmission power spectrum. This is in contrast 
to an estimate of the same quantity made by Zuo (Zuo 1993) who makes a prediction based on extrapolations 
of the column density distribution and of the number of absorbing clouds per. unit redshift (Zuo & Phinney 
1993). Zuo also assumes that the clouds are Poisson distributed, while our measurement incorporates the 
clustering in the IGM via our numerical simulation. 

An immediate problem presents itself: sightlines from which the mean transmission is measured are 
typically longer than the usual simulation box. We tackle this problem by expressing the variance of mean 
transmission in terms of the transmission power spectrum, and making use of a reasonable assumption about 
the behavior of the power spectrum on large scales. 

The mean transmission from one sightline is estimated using 

1 ^ 

^-mE^^ (6) 

1=1 

where N is the number of pixels, Fi is the observed transmission at pixel i, Fi = fi + rii where / = e~'^ is 
the true transmission and n is the noise fluctuation. We use the symbol F to represent the estimator, and 
/ to denote the true mean transmission. The variance of the estimated mean transmission is then 

4 ^ {F') {Ff = ^^[(F.F,) - {F){F,)] (7) 

hi 

where ct^ = (n^) is assumed roughly independent of position, and is the un-normalized two-point corre- 
lation of the transmission i.e. = (fifj) — f^, and Pf{k) is its one-dimensional Fourier transform. The 
symbol L denotes the comoving length of the spectrum from which the mean transmission is measured, and 
k is the comoving wavenumber. 

To evaluate ctt, we need to know the transmission power spectrum on scales generally larger than 
the size of the typical simulation box. It is expected that the transmission power spectrum takes the 
shape (not the normalization) of the linear mass power spectrum on large scales (i.e. essentially linear 



sin(A:L/2) 



kL/2 
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biasing; see Scherrer & Weinberg 1998, Croft et al. 1997, Hui 1999). We therefore use this to extrapolate 
the simulation Pf{k) to large scales (small fc's). We find that Pf{k) is well approximated by Pf{k) = 
B exp{~ak^) Jj^((ifc/27r)fcPmass(fc) where Pmass is the three-dimensional linear mass power spectrum. 

Becker et al. gave an estimate of ctt ^ 0.03 ± 0.01, / — 0.1, ai z = 5.5 using 4 different sightlines, each 
spanning Az = 0.2, which corresponds to L 57 Mpc/h. |^ An estimate of the noise term is provided by 
the error in the mean transmission, {o^/NY'-^ ~ 0.003. For the A^(fc = 0.03s/km, z ~ 2.72) = 0.74 case, 
the fitting parameters are B — 0.033 and a = 0.013 Mpc^/h^. Using eq. (||), we then find gt = 0.030 for 
A2(fc = 0.03s/km,z = 2.72) = 0.74, or = 0.031 for A^{k = 0.03s/km,z = 2.72) = 0.94, and (Jt = 0.028 for 
A^(fc = 0.03s/km, z = 2.72) = 0.58. The variance is similar between the different normalizations because 
each normalization requires a different in cq. (|l|) to match the mean transmission. This difference in 
Aa probably compensates for the effect of the different normalizations on gt- The predicted scatter of 
aT ~ 0.030 is consistent with the measured ctt of 0.03 ± 0.01. 

We apply the same methodology as the above to estimate ax z ^ 6. In Fig. ^ we show the results for 
a range of different Jhi's for each of our canonical power spectrum normalizations, (A^(fc = 0.03 s/km, z = 
2.72) = 0.58,0.74 and 0.94). ^ Photon noise is not included in the estimates of this figure. Even for 
relatively large Jhi's the scatter is small. For instance, for Jhi — 4.5 x 10~^, ctt = 1-1 x 10^^, assuming 
our fiducial normalization. This Jhi is large in that it already gives a mean transmission, / — 1.75 x 10^^, 
in excess of the observations. By Jhi = 1.4 x 10~^, the scatter is only ar — 2.8 x 10^^ for our fiducial 
normalization. The scatter depends somewhat on normalization, as one can see in the figure. To measure 
the scatter well would require data that are less noisy than the one discussed here, which has photon noise 
of (tT^j/iV)°-^ ~ 0.003, comparable to the predicted scatter. 

On the other hand, the smallness of this scatter makes it a possibly interesting diagnostic. As we have 
emphasized before, this predicted scatter ignores fluctuations in the ionizing background. For sufficiently 
small Jhi's, the IGM should be close to the epoch of reionization, and one would expect large sightline by 
sightline variations. An observed scatter well in excess of what is predicted would be an interesting signature. 



5. Discussion 

Our findings are summarized as follows. 

• The most stringent (1 a) lower limit on the neutral hydrogen fraction Xhi (eq. [§|) or upper limit on 
the ionizing background Jhi (eq. Q|) at z 6 is obtained from the observed mean Ly/3 transmission: 
Xhi > 4.7 X lO^"'. A comparison of this limit versus constraints at lower redshifts is presented in Fig. 
|l|. The fact that the neutral fraction increases by a factor of ~ 10 from rcdshift of 5.7 to 6 even though 
it changes by no more than a factor of about 2 from z = 4.5 to z = 5.7 suggests that z ~ 6 might be 
very close to the epoch of reionization. We emphasize that current constraints are still consistent with 
a highly ionized IGM at z ~ 6 - it is the steep rise in X^a that is suggestive of dramatic changes around 



^^The error on cry is estimated assuming Gaussian statistics and that the four lines of sight are independent. Then var{ax) = 
u?^/2n (see e.g. Kendall & Stuart 1958). 

^^The comparison across normalizations is done here at fixed Jhi while at 2 = 5.5 we compared the results of different 
normalizations at fixed mean transmission. We find that the dependence on normalization is larger at fixed Jjji than at fixed 
mean transmission. 
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or just before that redshift. We should also mention that the constraints on Xm are less subject to 
uncertainties in the IGM temperature compared to those on Jhi (see ^). 

• The existence of a long Gunn-Peterson {Lya or Ly/3) trough at z ~ 6, where little or no flux is detected, 
can also be used to obtain constraints on Xm or Jhi. This we call the fluctuation method: the fact that 
a long stretch of the spectrum exhibits no large upward fluctuations in transmission provides interesting 
information on the neutral fraction or ionizing background. The constraints obtained this way turn 
out to be fairly similar to those obtained using the mean transmission. We estimate that a reduction 
in noise by a factor of 4, or an increase in number of sightlines to 5, would result in constraints that 
are 2 times stronger (§^). 

• We develop a method to predict the dispersion in mean transmission measured from sightlines that 
are longer than the typical simulation box (eq. [Q and Fig. Our predicted dispersion is consistent 
with that observed at z = 5.5 (Becker et al. 2001). We also predict the scatter at redshift z — 6, which 
can be measured when more sightlines become available. Assuming a spatially homogeneous ionizing 
background, we predict a small scatter a,t z — 6, ax ^ a. few xlO~^, neglecting photon noise. The 
dispersion provides a useful diagnostic of fluctuations in the ionizing background - close to the epoch 
of reionization, one expects large fluctuations from one line of sight to another depending on whether 
it goes through regions of the IGM where percolation of HII regions has occurred. 

There are at least three issues that will be worth exploring. First, with more quasars at z 6 or 
higher discovered in the future, applying some of the ideas mentioned above would be extremely interesting, 
such as the measurement of the line of sight scatter in mean transmission, or the use of the Gunn-Peterson 
trough to obtain stronger constraints on the neutral fraction. Second, as we have commented on before, 
fluctuations in the ionizing background are expected to be important as we near the epoch of reionization. 
We have not discussed it here, but a calculation of the size of these fluctuations would be very interesting. 
Such a calculation will depend both on the mean free path of the ionizing photons as well as the spatial 
distribution of ionizing sources. The latter is probably quite uncertain, but useful estimates might be made 
(e.g. Razoumov et al. 2001). Lastly, a main source of systematic error which we have not discussed is the 
continuum placement. The mean transmissions at various redshifts given by Becker et al. are all obtained 
by extrapolating the continuum from the red side of Lya by assuming a power law of . The continuum 
likely fluctuates from one quasar to another, and therefore, it would be very useful to apply exactly the 
same procedure to quasars at lower redshifts where the continuum on the blue side can be more reliably 
reconstructed. This will tell us how much scatter (and possibly systematic bias) the continuum placement 
procedure introduces to the measured mean transmission. This kind of error is especially important to 
quantify given the limited number of quasars available for high redshift measurements at the moment. 
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z 


(/) 


4.5 


0.25 


5.2 


0.09 ±0.02 


5.5 


0.097 ±0.002 


5.7 


0.070 ±0.003 


6.05 


0.004 ±0.003 



Table 1: A summary of the observed mean transmission. The observation at redshift 4.5 is from Songaila 
et al. (1999). For this observation no error bars were provided by the authors. The observation at 5.2 is 
from Fan et al. (2000). The other observations arc from Becker ct al. (2001). Becker et al. (2001) have two 
observations at z = 5.5. The above mean transmission at z = 5.5 is the average of these two observations. 
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Fig. 1. — The top panel shows the neutral fraction of hydrogen at mean density as a function of redshift 
implied by the measurements of the mean transmission in the Lya forest. The point with the error bar 
pointing towards a completely neutral IGM comes from matching the mean transmission in Ly(3. The error 
bars include the 1 a uncertainty in power spectrum normalization and the 1 a error in the observed mean 
transmission. The dotted line is offered as a guide to the eye. It shows Xhi = 3.5 x 10~^(1 + The 
bottom panel shows the corresponding evolution in the ionizing background. 
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Fig. 2. — The upper panel shows the size of the jump in neutral fraction, Xm{z = 6.05)/Xhi(2: = 5.7), as a 
function of power spectrum amplitude. The amplitude is described by the value of A^(fc) = A-Kk^P{k)/{2'KY 
at 2: = 2.72 and velocity scale k = 0.03( km/s)~^. X]ii{z = 6.05) corresponds to the lower limit arising 
from the 1 a error in the mean transmission. The error bars come from the 1 a uncertainty in the mean 
transmission at z = 5.7. The lower panel shows the neutral fraction itself at 2; = 6.05. The dotted line is 
Xui = 5.8 X 10~''(A^(A;)/0.86)^-^, demonstrating how the neutral fraction scales with amplitude. 
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Fig. 3. — The upper panel a shows the probability Pa = J^p P{Fi...FN)dFi...dFN where Fi is the 
hya transmission at each pixel i of width aA. Here, N = 60, the noise per pixel is y/ip?) = 0.02 and 

Fm = 3-\/ (n^). The lower panel b shows an analogous probability Pp except that Fi now contains both Lya 
and Ly/3 absorption. Here, N = 48, ^ (n^) = 0.02 and 0.005 for dotted and dashed line respectively. The 
arrows indicate the corresponding 1 a upper limit on Jhi for these different probability distributions. 
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Fig. 4. — The lower panel shows the one-pixel (4A) probability distribution function of the true noiseless 
transmission (i.e. P(e~'^)de~^ gives the probability) for 3 different values of Jm- 0.004 (solid), 0.012 

(dotted) and 0.028 (dashed). The upper panel shows the probability distribution function of the noisy 
observed transmission F for the same three values of Jm- The negative values for F occur because of sky 
subtraction. 
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Fig. 5. — A prediction of the variance of the mean transmission, ctt, (see Section y) at z ~ 6, for several 
values of the ionizing background, Jhi- The estimate ignores contributions from photon noise. The triangles 
are for a model with power spectrum normalization A^(fc = 0.03s/km, z = 2.72) = 0.94, the squares 
A2(fc = 0.03 s/km, z = 2.72) = 0.74, and the pentagons A'^{k = 0.03 s/km, z ^ 2.72) = 0.58. 



